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ABSTRACT 


The need for direct turbulence measurements in the upper 
ocean arise from widespread requirements to correctly 
parameterize momentum and scalar fluxes across the air/ocean 
interface. Until recently these observations were limited by 
a lack of instrumentation capable of measuring the fine- 
structure velocity field down to dissipation scales. 

The recently developed CDV package allows simultaneous 
sub-centimeter resolution measurements of temperature, 
conductivity, pressure, shear and the 3 component velocity 
field. As part of the development of this system a Monte-Carlo 
simulation was used to analyze the performance of several 
spectral estimators of the mean acoustic Doppler shifts, from 
which the velocity components are derived. The selected 
algorithms were implemented on a Digital Signal Processor 
allowing real-time estimation of the velocity, shear and 
scalar quantities. 

To verify the performance of the CDV package,it was 
deployed off Wharf 2 in Monterey for 24 hours while 
simultaneously recording meteorological data. ’.<eynold's 
stresses, buoyancy fluxes and fine scale stratification were 
characterized and the surface gravity waves fiald identified. 
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I. 


INTRODUCTION 


Turbulence is a feature of flows and is characterized by 
its irregularity or randomness, diffusivity, large Reynolds 
numbers (highly nonlinear) and by three-dimensional vorticity 
fluctuations. Turbulent flows are also dissipative, requir’ng 
a continuous supply of energy, and limiting the smallest 
scales to be far larger than any molecular length scale. The 
smallest scale occurring in turbulent flows is called the 
Kolmogorov length scale [Ref.l]. These characteristics make 
any deterministic approach to the description of turbulent 
flow impossible, and observations very difficult. On the other 
hand, statistical approaches to the equations of motion lead 
to an unresolved system, requiring extra assumptions in order 
to obtain closure [Ref.z]. 

In order to characterize turbulence, we need to perform 
measurements of mean quantities and of fast varying 
quantities, requiring high spatial resolution, fast response 
sensors and stable observation platforms. A classical approach 
to separating mean quantities from turbulent fluctuation, is 
to assume a spectral gap, so that by using a running mean we 
are able to isolate perturbations from the mean flow (Reynolds 
averaging) [Ref.2]. 

A more complex procedure is required if one wants to 
separate a linear wave field from superimposed turbulence 


1 




[Ref.3 and Ref.4]. One method is just to subtract, from the 
fast varying variables, the components that are coherent to 
the wave field. The cross spectral variables and its coherency 
are used to identify the wave components and the correspondant 
frequency ranges. This method will only demodulate 
superimposed turbulence from the wave field assuming that no 
energy is being transferred (linear approach), and non-linear 
interaction between the two fields are not identified. 

Turbulence velocity observations are frequently limited by 
a lack of instrumentation capable of measuring the fine- 
structure temperature, salinity and velocity field, down to 
dissipation scales. The CDV package recently developed at the 
Naval Postgraduate School allows simultaneous sub-centimeter 
resolution measurements of temperature, conductivity, 
pressure, shear and the 3 component velocity field. These 
concurrent high resolution ocean microstructure measurements 
will aid our understanding of a range of phenomena that 
requires the determination of momentum or scalar fluxes. 
Unresolved problems include layer stability characterizations, 
turbulent closure problems, dissipation and energy cascade 
estimations, observation of entrainment and turbulence 
development, air-sea interaction, and the superposition and 
coupling of turbulence to mean advection and/or wave fields. 

As part of this work, real time signal processing 
programs, programed in C environment are developed, for the 
CDV package. This work is focused in creating and implementing 
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algorithms that estimate mean doppler shifts, from which the 
velocity components are derived; on the real time processing 
and filtering of the analog variables, including temperature, 
conductivity and shear; and on the continuous real time 
graphical monitoring of selected variables. Also, using a 
Monte-Carlo simulation, the performance of several spectral 
estimators of the mean acoustic doppler shifts are analyzed to 
guide the real-time processing algorithm choice. 

The CDV package, after being tested and calibrated, was 
deployed on the 2 and 3 August 1991, in shallow water off a 
Wharf in Monterey to test its ability to measure momentum and 
scalar fluxes. It was fixed to a 2 meter high bottom mounted 
tripod, with the instrument facing the surface at depths of 
the order of 0,5 meters, while simultaneously recording 
meteorological data. Small swell, weak winds and overcast sky 
were observed throughout the observation period. 

A procedure is developed that allows the estimation of 
Reynolds stresses and buoyancy fluxes based on the measured 
quantities. Spectral methods are used to characterize which 
part of the observed scalar variance fields are turbulent v.s. 
advected by the surface gravity wave field. 

This thesis is organized as follows. Chapter II summarizes 
the theorical background and motivations for microstructure 
observations. Chapter III describes a series of Monte-Carlo 
simulations that study the performance of several spectral 
estimators for the mean acoustic doppler shifts used to 
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estimate velocities. Chapter IV describes some recent methods 
and equipment that led to the development of the CDV package; 
the CDV package; the signal processing developed algorithms; 
and some validation tests and error analysis performed with 
real data. Chapter V describes the Monterey Wharf #2 
experiment presenting the Reynolds stresses and ocean fluxes 
estimation procedures. Finally chapter VI presents some 
conclusions and recommendations for future work. 
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II. MOTIVATION FOR MICROSTRUCTURE OBSERVATIONS 


Microstructure measurements are required for oceanography 
problems such as the study of mixing parameterizations, 
dissipation and energy cascade, surface layer entrainment, 
air-sea interaction and coupling, upwelling systems and mixing 
layer models. This chapter will introduce some theorical 
background to ocean turbulence emphasizing the experimental 
needs. 

A. REYNOLDS EQUATIONS 

The dynamics and thermodynamics of the ocean are described 
by a set of non-linear basic governing equations containing 
time and space derivatives requiring the definition of initial 
and boundary conditions [Ref.5]. Usually we do not have enough 
information to solve these conditions at smaller scales, so we 
must choose a cut-off eddy size below which we include only 
the statistical effects of these smaller scales. The resulting 
non-linear equations will be underdetermined with no 
analytical solution. Therefore, we are limited to 
approximations or numerical solutions after closing the system 
with additional assumptions [Ref.2]. 
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1. Basic Governing Equations 

The governing equations are the set of equations 
describing the general physical laws of fluid dynamics and 
their application to the particular case of ocean dynamics 
[Ref.6 and 2]: 

a. Equation of State 

This equation describes the thermodynamic state of 
water relating the several basic variables of density, 
pressure, temperature and salinity. A complete empirical 
relation in polynomial form is given, for example, by Millero 
et.al. 1980 [Ref.6]. However, this equation is difficult to 
handle and interpret. An approximation widely used is based on 
the fact that the water is nearly incompressible near the 
surface so that density changes of a constant volume are only 
due to temperature or salinity changes (Brian and Cox 1972) 
[Ref.5], 

p = p,-x [l-ax(r-T 3 ) +px(5-5,,) ] 

where a = 2139xlO''’K'^ is the thermal expansion coefficient and 
B = 0.751x10'^ is the salinity expansion coefficient and po = 
1025.96 kg/m’ is the density correspondent to the temperature 
To = 15°C and salinity So = 35 psu [Ref.5]. 
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b. Conservation of Mass (continuity) 


The conservation of mass for fluids, also called 
the continuity equation, states that the local rate of change 
of density is equal to minus the mass divergence, or: 


dt dxj 


ipUj) -0 


If the fluctuations of density are much smaller than the mean 
density it can be shown that the incompressibility assumption 
is valid and this equation becomes: 


dp _ ^^^7 

dt dx^ 


= 0 


i.e. non-divergent flow and constant density following the 
motion. 

c. Conservation of Momentum 

This balance describes Newton's second law of 
motion. If we assume that the vertical and horizontal scales 
are much smaller then the earth radius, this conservation law 
can be written using the tangent plane approximation [Ref.2]: 


dU. 

~dt 



-5,,g-€ 




dp ^ 1 

dx^ p dx^ 


1 






where: 


that: 


(5 = [1 Tn = n €=ri innq=123 231 312 

-1 m * n 1-1 innq=321 213 132 

0 otherwise 

Water can be assumed to be a newtonian fluid so 




"■ . 


8Xj 


where /i=l.4x10-3 kgm'^s"^ and are the molecular viscosity and 
the bulk viscosity respectively [Ref.2]. 

However if we assume that the fluid is incompressible and 
defining kinematic viscosity (v) as /j/p=l. 4x10*® m^s'^ we 

obtain: 

1 8x,. 

-_i_ =v-i 

p dx- dx/ 


If we assume that the motion is shallow, we can 
neglect density changes associated with the mass of the fluid, 
though we must retain their effect on the fluid weight, and we 
can further simplify the equations of motion in the so called 
Boussinesq approximation [Ref.2]: 


dU, 

it 


dU, 


P 


-c f n -1 


d^U: 

+ V-i 

8x,- 
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These equations represent the balance of the inertia and 
advection (LHS) by gravity, coriolis acceleration, pressure 
gradient and viscous stresses. 

d. Conservation of Heat 

The conservation of heat describes the first law of 
thermodynamics. If we neglect the change of internal energy 
due to viscous dissipation and assume incompressibility, it 
can be represented by: 


3-'T 

'-fc 


+ 


0 


r 


where the left hand side represents the storage and advection 
of heat, Kt=lxl0-7 m^s‘V(pCp) is the coefficient of thermal 
diffusivity, Q represents the rate of heat addition per unit 
mass by internal heat sources and Cp=3990 J.Kg’^K'* is the 
specific heat at atmospheric pressure and IS^’C [Ref. 5]. 

The term Q can be summarized by [Ref.2,5 and 6]: 

(1 -a) 


where: 


- short wave solar radiation (a = surface albedo) 
- net heat gain by long wave radiation 
Qp - net heat loss by long wave radiation 
Qp - net heat loss/gain by conduction 
Q.. - net heat loss/gain by evaporation or condensation 
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Q^. - net heat loss/gain by advection 
The terms Q* and will only be significant at the surface 
and can be estimated by the aerodynamical bulk formulas, based 
on bulk meteorological measurements [Ref.2]. 

The term representing the heat associated with the mass 
transport is usually parameterized by an eddy diffusion 
coefficient of the form: 


0,=K 


> d-T 

" 


The terms associated with the downward long wave 
and short wave radiation are more difficult to analyze. We 
must use empirical relations describing the penetration and 
absorption of radiation with depth as a function of 
wavelength. 

Finally the heat loss by long wave radiation can be 
estimated using Steffans law where (Qt « T*). 
e. Conservation of Salt 

Because of the form of the equation of state used, 
a salinity budget must be introduced of the form: 


as ,, ds d^S 


*S 


S 
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stating that the storage and advection of salt (LHS) must be 
balanced by the molecular diffusivity of salt and by the net 
salt gain/loss by internal sources. 

The term corresponds to the balance due to 
evaporation, precipitation and runoff effects. 

2. Reynolds Averaging 

The randomness found in the smaller scales of motion 
in a turbulent fluid makes a deterministic approach to solving 
these equations very difficult. If we assume that the energy 
at these scales is the superposition of a variety of small 
size, band limited eddies, we can separate these components 
from the other mean effects by averaging over a sufficiently 
large interval corresponding to a stationary low frequency 
energy spectral component. This concept presumes the existence 
of a spectral gap in the component variance spectra [Ref.2]. 
However, the presence of internal waves, surface gravity waves 
and many other phenomena, does not always allow us to assume 
the existence of this gap between forcing and dissipation 
frequency ranges. Nevertheless, we can assume that each 
variable is composed by a mean and a perturbed part, which 
averages out to zero. However, the turbulent part of the 
variance will have non-zero covariances with the other 
turbulent components. 


11 







starting from the previous basic governing eguations 
and after separating the mean and perturbed components and 
averaging, the resultant set of equations is: 

1. equation of state 

p = P3x[l-ax(T-ro) ^pxis-s,)] 

2. conservation of mass (continuity) 

dVZ _ du'- 
dx- dx. 


3. conservation of momentum 


du. _ du. 

=-5 , .a-e . 


dt 


dx- 


k ^ ck 


7r_4 JP 


d'-^lL 


dx. dx. 


dx. 




4. conservation of heat 


dt ' dx, ' dxU 


U' 


5. conservation of salt 


at ^dx,^ ^dx? 



12 






The new terms introduced by this procedure are the 
divergences of momentum, temperature and salinity fluxes, and 
they represent the interaction between the mean and the 
perturbed flow. For these terms to be significant, we must 
have a large Reynolds number ( = advection terms / 
molecular viscosity terms), so that it is usual to neglect the 
viscosity terms in this set of equations. The set of equations 
represent an unresolved system of 7 non-linear equations with 
22 unknowns and are known as the Reynolds equations. 

Before looking for any type of solutions we first must 
close the system by finding prognostic equations for the new 
variables in the form of fluxes and variances, eventually 
reducing the number of unknowns. To obtain the forecast 
equations for the fluxes and variances we must first find 
prognostic equations for the perturbatjons. Following the 
procedure described in Stull, 1988 [Ref.2], this can be 
achieved by subtracting the Reynolds equations from the 
governing equations, resulting in the system: 


dll': _ dll': U, , du', 




dx. 


dx. 


P ' P 


1 dp ' dru'.^ Q — 1 -J- 

--i + - {U .U,) 

dx/ dx, ^ ' 
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bt ^ ax. 


+ u 


bT ^ / bT' d^T' 

-T - +a --:r-A>- 

bx. - bx- dx/ 


c„ ax- 


bS\jr dS' 
bt - bx. 


■u 


bs 

ax. 


■u 


bs' 

ax. 


=i<r. 


b'^s' 

bx,^ 




bx. 


( 17 ^) 


In order to obtain equations for the reynolds stresses and for 
the scalar fluxes, the following procedure was developed by 
Stull, 1988 [Ref.2]: 

1. Multiply the equations for a generic variable q', by u'^, 
and Reynolds' average. 

2. Sum these resultant equations, obtained by Reynolds' 
averaging the perturbed momentum equations multiplied by q,'. 

3. Use the continuity balance to transform advective terms 
into flux form. 

4. Neglect the coriolis terms and the pressure, radiation and 
molecular diffusion. 

The resultant equations, are; 




-r ,.u-. 

bt ^ dx, - dx, ^ ^ dx,- dx 
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dY'u ', _ dT’i:'. 


dt 




dx 


^ dx 


-u ,u'. 


dx. dx. 


(i u',u ^) 


■6,,(W>9^W^-2e^.r 

p P OX. 


ds'i:'. _ ds'u^ 

J.+U, 


dt 


-S 


T-^^i 




ds 


ax, ax, ax, ax^ 


iS'u’jU^,) 


P 


dx. 


u ,S 


As we can see, these equations do not close the system because 
they introduce new variables in the form of statistical third 
moments. If we proceed looking for prognostic equations for 
these third moments, they will introduce in a similar way 
fourth order type moments and so forth. 

Therefore in order to be able to use these Reynolds 
equations additional assumptions must be used. 

3. Turbulent Closure 

If there was no turbulence and the statistical second 
moments of the Reynolds equations were negligible, they would 
form a closed system of equations describing the mean motion. 
This is typical of the deep ocean, where turbulence events are 
episodic and of a low mean level [Ref.6]. 
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However, these second moments are generally not 
negligible on the upper layer of the ocean, where atmospheric 
fluxes can be strong. The upper layer, where the Reynolds 
stresses and other turbulent fluxes are important, is known as 
the turbulent ocean mixing layer. It is not necessarily equal 
to the well mixed upper layer of nearly constant density known 
as the remnant mixed layer [Ref.7]. Therefore, in this study, 
our closure problem is essentially a boundary layer problem 
where mixing processes and air-sea interactions take place. 

In order to close the system of equations, we must 
first decide to what order we want to go to in our prognostic 
equations, and at that point parameterize the statistical 
moments based on known variables. Hence parameterizations are 
an important tool that rely essentially on observations of 
slow varying variables. However, a complete study of these 
parameterizations, to see how the observed variables relate to 
the fast varying ones and to validate the procedure, will 
involve the simultaneous observation at a wide range of 
scales. Until recently these type of observations suffered 
from a lack of instrumentation capable of providing a complete 
picture of the ocean microstructure, similar to those that 
exist for micrometeorology observations. 

The turbulent closure approximation may be made to 
order 0, 1, 2, etc, depending on the highest order moment in 
the prognostic equations. Some theories only use a portion of 
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the available equations introducing what is known as the half¬ 
order closures. 

Two major approaches of turbulent closures are used 
that are applicable at any order [Ref.2]. In local closure, we 
parameterize unknown quantities at one point by values or 
gradients of known quantities at the same point. It assumes 
that turbulence has a similar nature to molecular diffusion. 
This is like K-theory, where we assume the stresses are 
directly proportional to the mean shear. In the other 
approach, known as nonlocal closure, an unknown quantity at a 
particular point is parameterized by known quantities at many 
points in space. It assumes turbulence as a superposition of 
eddies transporting fluid and its properties, like an 
advective process. Examples of nonlocal first order closure 
like the Transilient Turbulence Theory and the Spectral 
Diffusivity Theory are described in Stull, 1988 [Ref.2]. 
Usually higher order local closures and nonlocal closures 
yield more accurate results, but they do so by adding 
computational expense and complexity [Ref.2]. 

1. Local Zero-Order Closure 

In this approach, no prognostic equations are 
retained, not even equations for mean quantities. Mean flow, 
temperature and salinity are parameterized directly as 
functions of space and time. Adimensional variables are 
usually defined for certain classes of problems. 
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2. Local Half-Order Closure 


This closure scheme uses only a subset of the first 
moment equations. A variation of this method yield the so 
called Bulk Models, where temperature profiles are assumed to 
have a certain shape, exploring homogeneity in the mixing 
layer and using vertically integrated equations [Ref.8], 

3. Local First-Order Closure 

In this scheme we retain the prognostic equations 
for the mean quantities and we parameterize the statistical 
second moments (variances and cross-variances). 

For a generic variable ”1" one possible 
parameterization of its turbulent flux will be: 


u',l'=-Kj 


dL 

dx. 


where Ki with units m^s'^is known as eddy viscosity for the 
velocity case and as eddy diffusivity coefficient for the 
scalar variables. This approximation is called K-Theory and 
frequently fails whenever large eddies are present or when 
there is not a dominant scale [Ref.2]. A useful interpretation 
of these approaches can be made by looking at analogies to the 
viscous stress expressions for Newtonian fluids. We can 
express the Reynolds stresses in terms of the shear, by 
substituting the molecular viscosity coefficient with an eddy 
viscosity coefficient, obtaining: 
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dU 


^reyno2ds P^in 


Since turbulent mixing is much more effective than 
viscosity, the eddy viscosity should be much greater then the 
molecular viscosity [Ref.2]. Also while molecular viscosity is 
a property of the fluid, the eddy viscosity is a property of 
the flow so that this coefficient should vary as turbulence 
level varies. Good parameterizations should rely on the Monin- 
Obukhov length scale, on the Richardson number (Ri) and on the 
stratification, as we can see in the examples presented in 
Stull, 1988 [Ref.2]. 

Another useful interpretation of this theory is the 
so called Mixing-Length Theory derived by Prandtl (1925) 
[Ref.2] for neutral conditions. If a turbulent eddy moves a 
water parcel upward by a distance z' and if we assume that 
there is no other mixing or changes on "1" within the parcel, 
then the perturbation of “1" can be written as: 


7 /= 7 -7 =- ( ^ ) 

suzrou.idir.<J -^parcel ' ' 


similarly 


/ 

^ ~^suizoundipg ^pazcel~ ^ ^ ^ 
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In order to have upward motion we must have some vertical 
velocity w'. If turbulence is such that it can be assumed that 




c. u 


w^=-c. u' 


for 

for 


dz 

ffu 

dz 


<0 

>0 


where "c" represents a positive constant. Then we can define 
mixing length "Ml" as 

Ml ~-cz'^ 


A vertical eddy flux of the scalar variable "1" will have the 
form: 


cz' dz 


so that 

dz 


Therefore Ki should increase as shear increases or turbulence 
increases and as the ability of turbulence to cause mixing 
increases. Finally a remark should be made about the 
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assumption of linear gradients valid only over small 
distances, noting that this theory is a small-eddy theory. 

4. Turbulence Kinetic Energy (TKE) 

TKE measures the intensity of turbulence and is 
directly related to the momentum, heat and salt transport 
through the mixing layer. Its balance is important in the 
physical interpretation of the fluxes that are present on the 
Reynolds equations, and in describing the physical processes 
that generate and maintain turbulence. The total kinetic 
energy per unit mass (KE) is given by: 



using Reynolds averaging we can identify the following 


components on the mean KE: 


K=Y' —1''" KE of Che wean flow 
; =1 


— u'l KE of turbulence {TKE) 

2 


From the Reynolds equations we can obtain the budget for the 
KE of the mean flow on the form [Ref.2]: 
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^.TT. 

dt 


dK - -rr d 


= -p.£7 


dx. 


{u',u‘ 


■-AI7-|P 

p 


_ _ a 

-Wg-^L\\ 


U. 


dxt 


Interesting features in this equation are the fact that 
Coriolis has no effect on the TKE, showing that it is only a 
redistribution term. Another important feature is the explicit 
form of the Reynolds stresses representing the transformation 
of energy between the large scales of the mean flow and the 
smaller scales of turbulence. 

From the definition of potential energy (tp) and by 
using the Reynolds averaging procedure we can obtain the 
following potential energy budget: 


dt 


=Wq->-^w'p 'q 
P 


We can see from this budget that the scalar turbulent fluxes 
represent the transformation between the mean potential energy 
and turbulence. 

After identifying the role of the Reynolds stresses 
and scalar fluxes it is important to see how they balance, 
interact and develop in time. This can be accomplished by 
establishing a TKE budget. Using the prognostic equations for 
the perturbations and using an analogous procedure to that 
wich we used to obtain the prognostic equations for the 
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Reynolds stresses and scalar fluxes, we can obtain prognostic 
equations for the variances of the perturbed velocity 
components that yield the following TKE balance [Ref.2]: 


dt 



p ^ ^ 




dx. 


d(u',k) 


1 d (u'vp') 

■p 8 a'_- 


The LHS of the budget represents the storage and advection of 
TKE by the mean flow. This latter term is often assumed 
negligible by assuming homogeneous turbulence over the mixing 
layer. 

The first term on the RHS represents the buoyancy 
production/consumption. Corresponds to the transformation of 
TKE into potential energy and vice versa. Usually it is 
positive during the day and negative during the night due to 
diurnal solar heating. 

The second term represents the mechanical or shear 
production and is directly related to the Reynolds stresses, 
reflecting any energy sources due to stirring or mixing and 
the interaction with the kinetic energy of the mean flow. 

The third term represents the eddy transport of TKE, 
acting as a source or loss whenever there is a convergence or 
divergence. 

The fourth term represents the pressure correlation 
term, describing how TKE is redistributed by pressure 
perturbations. It is this term that is responsible for the 
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interaction of turbulence with internal or surface gravity 
waves. The small scale perturbations can be due either to 
waves or turbulence and from linear wave theory (u'^p') 
represents energy flux or power; therefore we can expect 
interaction between these two phenomena to exist by sourcing, 
sinking and redistributing TKE. 

The last term represents the viscous dissipation rate 
or the transformation of TKE into heat. 

B. MIXING LAYER ANALYSIS 

Mixed layer depth can be defined as the depth at which 
density increases by some predetermined amount above the 
surface value, while mixing layer depth can be defined as the 
depth at which the dissipation rate drops monotonically below 
some arbitrary small value [Ref.8]. This definition for mixed 
layer is useful for observational purposes, but do not reflect 
the depth range of mixing, i.e. the layer that is dynamically 
unstable. Static stability is a measure of the capacity for 
buoyant convection that is independent on the flow and based 
only on density profiles. Dynamic stability measures the 
tendency for mean flow to break into smaller scales, even 
under statically stable conditions [Ref.6]. While statically 
unstable flows become or remain turbulent, statically stable 
flows may be laminar or turbulent depending on many factors 
that can be interpreted based on the TKE balance. 
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We can view the dynamic stability in terms of a balance 
between the two different effects of stratification and shear 
whenever the inertia is much greater then the viscous effects 
(high R,) and the reynolds stresses and scalar fluxes are 
significant. While stratification will tend to damp out any 
perturbation, the shear of the mean flow will tend to enhance 
perturbations by supplying the required energy. 

A description of a mechanism for the dynamic generation of 
turbulence is summarized by Stull, 1988 [Ref.2] based on 
laboratory experiments. He describes the following sequence of 
events: 

1. shear across a density interface (laminar flow) 

2. a critical value of shear is reached and flow becomes 
dynamically unstable, generating gentle waves at the interface 

3. waves grow in amplitude and eventually break 

4. between waves crests lighter fluid entrains the denser 
one and patches of static instability results 

5. this static instability combines with continuous 
dynamical instability causing each wave to become turbulent 

6. turbulence at the interface causes mixing of the two 
different fluids (entrainment) 

7. mixing can reduce shear below critical value and 
eliminate dynamical instability 

8. if there is no continuous supply of energy to restore 
the shear, turbulence will decay and flow will become laminar 
again. 
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This conceptual model for generation of turbulence seems 
consistent, however due to the episodic and 3 dimensional 
nature of turbulent events in the ocean, observations of all 
the processes are very difficult to perform. Also the lack of 
microstructure instrumentation able to directly observe the 
velocity field is a major limitation on our understanding of 
these processes. 

Usually we rely on nondimensional parameters to compare 
these factors, like the Reynolds number and the Richardson 
numbers. The Reynolds number compares the inertial and the 
viscous terms in the equations of motion, evaluated when the 
Reynolds stresses are more important then viscous stresses. 
This happens only when the non-linear effects are significant. 
For turbulence to exist, it is therefore necessary to have a 
high Reynolds number, but this by itself it is not a 
sufficient condition, because it makes no requirements about 
the supply of energy necessary for turbulence to persist. 

The Richardson number compares the buoyancy production to 
the mechanical production terms of the TKE balance, and 
evaluates when the shear will be able to overcome 
stratification to allow dynamical instabilities to exist 
[Ref.2]. Usually we have high Reynolds numbers in the ocean; 
therefore, the Richardson number is the critical criteria to 
determine if the flow is turbulent or not. 

The flux Richardson number is defined as [Ref.2]: 
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A r-*- 

(u'.uV) 


dU^ 

dx. 


For Rf less then zero we have static instability (hence 
turbulent) flows. Rj equal to zero corresponds to neutral 
flows. The critical conditions below which flows are turbulent 
is R, = 1. 

Using K-theory, we can parameterize the fluxes based on 
the mean shear and obtain the gradient Richardson number 
[Ref.2]: 


R,= 


p dz 



dv 

dz 


)^] 


The critical value below which laminar flows become turbulent 
is accepted as 0.25, while the critical value for turbulent 
flows to become laminar is 1 [Ref.2]. This hysteresis is due 
to the fact that, for turbulence to start, we need the 
instability and a trigger mechanism; while instabilities will 
start at Rj near 1, the trigger mechanism will only occur near 
R, = 0.25 [Ref.2]. 
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Because usually we are limited to discrete observations at 
a certain depth, the bulk Richardson number is defined as 


[Ref.2]: 

^ p[Alr-Av^] 

This discretization will change the critical values but the 
same critical value of 0.25 is generally used. 

From this discussion we can conclude that direct 
microstructure observations are very important on the 
characterization of the mixing layer. Mixing processes have a 
special importance in coastal processes, controlling the rate 
of cooling or warming of the surface waters and the dispersion 
or diffusion of other scalar quantities. The role of vertical 
turbulent mixing on the coastal upwelling systems is 
determinant in the characterization and understanding of these 
type of phenomena, as we can see in the paper by Brink (1980) 
[Ref.9]. Therefore mixing layer observations, are need in 
order to study: mixing layer deepening or entrainment; and air 
sea interactions, as we can see in the papers by Dewey et.al. 
and Mourn et.al. [Ref.7 and 10], characterizing coastal eddies 
and the enhancement of ocean fronts by turbulent mixing, and 
in the paper by Uwe et.al., analyzing upwelling relaxations 
during the CODE experiment [Ref.11]. 
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C. SPECTRAL DESCRIPTION OF TURBULENCE 

The spectrum, defined us the Fourier transform of the 
autocorrelation function, gives a picture of how variance or 
energy is distributed over frequency or wavelength. An eddy of 
a certain wavenumber "k" can be defined as a linear 
superposition of the disturbances containing energy in the 
vicinity of "k" [Ref.l]. These eddies will be exposed to the 
strain-rate field of larger ones, and nonlinear interactions 
and consequent energy transfers are expected. As a resul'*', the 
eddies will loose their identity within one or two periods or 
wave"jngths. 

The Fourier Transform is a decomposition into wave 
components of different wavelengths associated with the 
fourier coefficients. An eddy, however, has many fourier 
coefficients changing in time. Therefore more sophisticated 
transforms and a nonstationary approach should be used if we 
want to decompose the velocity field into individual eddies. 

When analyzing the wavenumber spectra, i.e. how energy or 
variance is distributed in scale, we can expect to find 
contributions down to the sub-centimeter scales of 
dissipation, after which the mechanical energy is transformed 
into heat by viscous processes, and the energy drops down to 
zero [Ref.l]. If our data is not sampled below these scales 
the spectra will be aliaised. Therefore to understand both 
momentum fluxes and dissipation we must measure a wide range 
of scales, including sub-centimeter dissipation scales. 
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a. Velocity Spectrum 

Large eddies are affected by the mean flow, such 
that they have a steady anisotropy, whereas no permanent 
anisotropy is induced at an appropriately small scales. As the 
eddy size becomes smaller, we can expect that the strain-rate 
field becomes isotropic in a mean sense. This concept is 
called local isotropy and it can not exist if the Reynolds 
number is not large enough (at least 100 [Ref.l]). When local 
isotropy prevails, the time scales are much shorter then those 
of the mean flow, so that small eddies respond quickly to 
changing conditions. This area of the spectrum is called the 
equilibrium range and includes all high wavenumbers. In this 
area of the spectra the scales are too small to receive energy 
directly from the mean flow. Because energy is finally 
dissipated, the amount of energy arriving to this scale should 
correspond to the total energy being dissipated by viscosity. 
The smallest scale that can be reached on this range is called 
the Kolmogorov microscale defined as [Ref.l]: 



where v is the kinematic viscosity coefficient and e is the 
dissipation rate given by [Ref.l]: 
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e = 2.vfk^.S.dk 
0 


and "S” is the observed velocity spectrum. 

For small wavenumbers, viscosity is not important. 
The important parameters to be considered are those relating 
transfer of energy from the mean flow to turbulence, and the 
energy transfer to the smaller scales. The area of the 
spectrum corresponding to the superposition of theL,e scales 
with the Kolmogorov spectrum is often called the Inertial 
Subrange and exists only for large Reynolds numbers (bigger 
then 10®). At this range, no energy is added or dissipated, so 
that the energy flux across each wavenumber is constant. It is 
characterized by a slope of -5/3 on a loglog plot and is 
scaled by dissipation. Therefore, from observed spectral 
estimates in this range, we can estimate the total 
dissipation. 

Spectral estimates will depend only on the 
wavenumber and on the total dissipation. Using similarity 
theory [Ref.2], we can expect the spectral estimate to be: 

S(k) =a,j,.e^^Kk--'^ 

[5] =Tn^s 
where [e] 

[k] 
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where a„ is known as the Kolmogorov constant with ranges from 
1.53 to 1.68 [Ref.2, pp.390]. 

For time spectra of velocity, two cases should be 
considered. The Lagrangian time spectra related to the 
temporal evolution seen by an observer moving with the 
turbulent velocity fluctuations and the Eulerian spectrum 
related to the evolution at a fixed point where the mean 
velocity is zero. Their shape is usually different due to 
advection, but we can still define the inertial and 
equilibrium sub-ranges in a similar way. 
b. Temperature Spectrum 

The temperature spectrum is defined in a similar 
way to the velocity spectrum, and shows the same nonstationary 
character, with energy being added, dissipated and transferred 
between different scales. If the Reynolds number is high 
enough, there is also an equilibrium range, an inertial 
subrange and an inertial-convective subrange where the 
temperature fluctuations are simply convected and where the 
spectral transfer of temperature variance is constant. In 
certain cases when thermal diffusivity becomes important 
within the inertial sub-range and where viscosity is not yet 
important, there is a Inertial-Diffusive subrange; this occurs 
in fluids with small Prandtl number (ratio of kinematic 
viscosity to heat conductivity). For large Prandtl number 
fluids, like water, viscous effects occur first and then 
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thermal diffusive effects; therefore we can also define 
Viscous-Convective and Viscous-Diffusive subranges, 
c. Batchelor Spectrum 

The Batchelor spectrum corresponds to a theorical 
shape of the wavenumber spectrum of the temperature gradient 
fluctuations for homogeneous, locally isotropic turbulence in 
a homogeneous fluid with a large Prandtl number [Ref.12 and 
1]. It is defined by the temperature variance dissipation 
rate, molecular diffusivity and by the Batchelor wavenumber 
defined as: 

iCj, = ^ 


This wavenumber as been identified as the cutoff wavenumber 
[Ref.12] . 

The effect of stratification was studied by Dillon 
and Caldwell,1980 [Ref.12], who proposed that its effect 
should be more noticeable in the low viscous-convective 
wavenumber ranges when turbulence is weak compared with 
stratification. One of the most important applications of this 
theorical spectrum is the estimation of the kinetic energy 
dissipation rates from the spectral cutoff wavenumber, as we 
can find in the work of Dillon and Caldwell,1980 [Ref.12]. 
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III. DOPPLER SPECTRUM ESTIMATION 


The bistatic velocity estimation that will be discussed in 
the next chapter, requires a fast but accurate doppler 
spectral estimator. This estimator must rely on a small number 
of data points and be reliable even for weak signals submerged 
in strong noise. Several spectral estimators will be analyzed 
and Monte-Carlo simulations performed in order to study their 
performances. 

Spectral estimation corresponds to the estimation of a 
function that is, by definition, the Fourier Transform of the 
autocorrelation function. The performance of a spectral 
estimator should be based on how well we are estimating this 
function over a continuous range of frequencies [Ref.13]. 

The use of Spectral Estimators to detect the presence of 
signals superimposed in white random noise demands a trade off 
between detectability and confidence of the detected values. 
Usually confidence is obtained by dividing the whole data set 
into smaller subsets, estimating the spectrum of each subset 
and then averaging the ensemble of estimates. If we consider 
that each spectral estimate corresponds to 2 Fourier 
coefficients and therefore has 2 degrees of freedom, the 
ensemble average based on N subsets will increase the number 
of degrees of freedom to 2N, improving the viability of the 
overall estimate. However the spectrum based on these subsets 
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will have less resolution and the spectral estimate will be a 
smoothed but more accurate version of the unsegmented one 
[Ref.14]. 

There are applications where we need only the behavior of 
the spectral estimator at a specific frequency or in a 
specific region, such as the frequency estimation of a 
sinusoidal signal superimposed in white noise. In this case, 
the smoothed version of the spectrum will be a poor estimator 
of the sinusoidal frequency as we are seeking a parameter 
estimation (and not spectral estimation), that involves 
different statistics [Ref.13]. 

The frequency estimator based on a finite number of data 
points represents an estimate of the true frequency that we 
would have if we had an infinite number of points. Therefore 
we need to characterize, in statistical terms, the behavior of 
these estimates. In this work, the Bias and Variance of the 
frequency estimates will be the statistics that we will use to 
compare the performance of several estimators. 

Because determination of these statistics for a wide 
variety of algorithms is not analytically treatable, Monte 
Carlo simulations are used to obtain estimates of these 
statistics in terms of the Signal to Noise Ratio (SNR) of the 
input signal and of the record length used in the estimator. 
The results of this simulation will be used to determine the 
appropriate doppler frequency estimator, for the bistatic data 
stream. 
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A. ESTIMATION OF A SINGLE COMPONENT COMPLEX SIGNAL 
SUPERIMPOSED ON GAUSSIAN WHITE NOISE 
1. Peak FFT 

This approach can be interpreted as an estimation of 
the parameters of a single component sinusoidal signal of 
frequency f, superimposed in gaussian white noise, of the form 
[Ref.15]: 

x(/7) =A EXP{j2Tlfn + <^) *zin) 

where z(n) corresponds to complex gaussian white noise and 
and "f” (amplitude ,phase and frequency), assumed 
constant but unknown, are the parameters to be estimated. 

The probability distribution function (PDF) of x(n) 
corresponds to a shift in the multivariate gaussian 
distribution function associated with the white noise, and has 
the form: 


=—-- EXP[ - iX-£) ] 

n-Met 

where "R„" is the noise correlation matrix and "s" the 
sinusoidal signal. Since we are assuming white noise, "R„" 
will be a diagonal matrix. The sinusoidal signal can be 
represented by the vector: 

£=A EXP(j^) [1 EXPijZUf) EXPijZUfZ) . . . EXPijZUfiN-l)]'^ = 
EXP[j^) £ 
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where e is a vector containing the values of the complex 
exponential of frequency f, at the several values of n. 

To obtain a maximum likelihood estimator (MLE), we 
need to find the maximimum of the PDF, which will be 
equivalent to minimize the hermitian form: 

■S' = = (K~A EXPij^) ii^-A EXP(j<f?) s.) 

First, if we minimize in terms of the complex amplitude 
Ac=A EXP(j(p), we obtain [Ref. 15]: 

A. = ^ 


which corresponds exactly to the discrete fourier transform of 
the data evaluated at f (assumed known in this case). Now if 
we substitute this expression into the hermitian form, we 
obtain: 

N 

to minimize this expression in terms of the frequency we need 
to maximize: 
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N 




■V-J 


= — ly^ x(n) EXP{-j2'n.fn) P 


Therefore the MLE of the frequency will be the frequency where 
the square of the absolute value of the discrete fourier 
transform (periodogram) attains its maximum. 

In this work, we used a 128 point FFT algorithm to 
estimate the periodogram where the peak frequency will 
correspond to the MLE estimator. 

This procedure will limit our resolution to 7.8x10'^ 
of the sampling frequency, so an interpolating procedure 
should be used in order to increase the estimator resolution. 
Nevertheless this estimator suffers from the fact that, for a 
small number of points, discretization errors are highly 
emphasized and its performance degradated if one wants to 
obtain faster estimations. 

2. Pulse Pair 

Traditionally, the Pulse Paj.r algorithm is used to 
estimate the mean frequency of a time series with a single 
frequency component signal, which corresponds to the frequency 
of maximum energy [Ref.16]. It is based on a very simple and 
fast algorithm that has given accurate results even for data 
sets with a small number of points. However, white noise 
contamination will induce a deviation from the peak value 
toward the zero frequency as the SNR decreases. 
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The algorithm uses the direct evaluation of the phase 
from the autocorrelation function estimate. In this project 
only the phase of the first lag autocorrelation function was 
used and computed using a classical approach based on 
calculating covariances of two consecutive points averaged 
over the sequence length. Other valid correlation estimate 
methods could be used such as the Modified Covariance Method 
[Ref.15], which has the property of giving estimates without 
running out of the data set. It will be shown later, that this 
method is equivalent to an eigenvalue decomposition of an 
order 2 correlation matrix (MUSIC -order 2 fast algorithm). 

Higher lag autocorrelation estimates can be used 
either to improve the estimates when the frequencies are much 
smaller then the sampling frequency, or to obtain more 
credible estimates. However this procedure may create problems 
with the uniqueness of the estimate, and will be limited by a 
lag, corresponding to the Nyquist rate of the highest 
frequency present on the detected signal, after which aliaised 
frequency estimates occur [Ref.16]. 

3. Zero Cross 

This method consists of the evaluation of a dominant 
frequency by the number of zero crosses or the number of sign 
changes on a zero mean time series. When working with complex 
signals, we must start by defining criteria for a zero cross. 
A point in the real domain corresponds to a line parallel to 
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the imaginary axis in the complex plane. Then using this 
generalization, we can define a zero cross (real) as the 
change of sign on the real components of the time series. On 
other hand, a point in the imaginary domain is a line parallel 
to the real axis in the complex domain; then a zero cross 
(imaginary) corresponds to a change in sign in the imaginary 
components of the time series. Therefore, we can define a 
complex zero cross as a change in sign either in the real or 
imaginary part between consecutive points; therefore the zero 
crosses can either be determined in the real or imaginary 
parts. To get a better estimate we can also average both 
estimates. 

By itself, this simple method is highly sensitive to 
high frequency noise contaminations that will increase the 
number of zero crosses. This effect can be partially 
compensated by choosing a hysteresis after each cross in order 
to minimize the high frequency contaminations [Ref.17]. In 
order to apply this technique, we must know the maximum 
possible frequency that we can have in the signal after which 
no frequency will be detectable. 

Another approach developed by B.Kendem [Ref.18] 
suggests the incorporation of low-pass filters and periodogram 
analysis to obtain more robust estimates. This procedure was 
not developed in the present analysis, though some interesting 
results may be obtained. 
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4. AR Modeling (Burg's method) 

AutoRegressive (AR) modeling is one general class of 
parametric models that we can use to study time series when 
some previous information exists [Ref.15]. It is based on the 
assumption that the value of the record at a particular time 
is determined by a linear combination of P past values, where 
P corresponds to the model order, superimposed in white noise 
(u[n]), in the form: 

P 

x[n\ =-^2 ^ x[n-k] +u[n] 

k-i 


The model coefficients are usually determined based on the 
autocorrelation function of the series using any of the 
estimators (usually the modified covariance method gives the 
best results) by solving the YULE-WALKER equations. 

Another method to compute the model coefficients is 
based on an extrapolation of the autocorrelation function that 
has maximum entropy, so that the effects of truncation 
(spectral smearing) are minimized [Ref.15]. The estimated 
coefficients are then computed based on this assumption 
directly using Lagrange Multipliers. This method is known as 
BURG'S technique and is the one that was used in the present 
study [Ref.15]. 

The spectral estimator will correspond to the z- 
transform of x[n] evaluated on the unit circle, or; 
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where is the sampling frequency, H(z) the z-transform of 
the signal x[n], A(f) the denominator polinomial of H(z) after 
the transformation z=EXP(-j .2.7r.f), a the model coefficients 
in vector form, and p„ the noise power. 

The spectral estimators will have singular points for 
a discrete type spectra (pure sin'isoids) corresponding to 
zeros on the denominator of H(z); therefore, in our case of 
only one pure complex sinusoid in white noise, we can just 
look for the frequency that will be closer to a singularity. 
This can be accomplished by looking for the zeros of the 
denominator of the general spectral estimation expression and 
choosing the one that is nearest to the unit circle. In this 
form, the frequency estimation will be totally determined by 
the AR model coefficients. This is the method we used for the 
AR dominant frequency estimation. Because of its complexity, 
we should expect this method to be computationally intensive, 
and the results will also depend on the root finder that we 
use. 

5. MiniiBuiD Variance Method 

This method introduced by Capon [Ref. 13] is also 
called the Maximum Likelihood Method. It will describe the 
distribution of variance over frequency, but the area under 
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this estimate will not corret >ond to the total variance 


[Ref.13]. Therefore it is not a true spectral estimator, but 
is particularly suitable for single frequency estimation. The 
algorithm computes the spectral estimates in accordance to the 
output power of a set of narrow band finite impulse response 
(FIR) filters whose bandwidth or resolution are dependent on 
the frequency in an adaptative manner (in comparison with the 
periodogram that has a fixed resolution). The filter 
coefficients are computed so that the output variance is 
minimum and its gain is 1 (no distortion). Ir, fact we do not 
need to compute all these coefficients explicitly because the 
spectral estimators can be derived directly from the inverse 
autocorrelation matrix in the form [Ref.13]: 


P 


MT 



_ 1 _ 


where R corresponds to the Toeplitz autocorrelation matrix, 
that can be estimated in several forms. In the present work 
the modified covariance method was used. 

It can be shown that it gives spectral estimates that, 
for a given order P, correspond to the average of the spectral 
estimates of the AR models of order 1 to P [Ref. 13], so its 
performance should result in a smaller variance. Its 
resolution is between that obtained by the classical spe tral 
estimators and by an equivalent AR model. Therefore, in order 
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to obtain the single frequency estimate, we just need to 
compute the maximum of the above expression. However, due to 
its complexity and dependence on the resolution used in the 
maximum value determination, we should expect this algorithm 
to be even more computational intensive than an equivalent 
order AR method. 

6. Multiple Signal Classification (MUSIC Algorithm) 

This algorithm is one of several that analyze the 
column space of the autocorrelation matrix into a signal 
subspace and an orthogonal noise subspace using a Singular 
Value Decomposition (SVD) [Ref.15]. This is accomplished by 
separating the matrix into its eigenvectors so that the 
eigenvectors with the N largest eigenvalues, span the signal 
subspace and the others the noise subspace, where N is the 
number of complex components of the signal. 

This estimate, like the minimum variance method, is 
not a true spectral estimator because it does not preserve the 
total variance, however it will work as a discrete frequency 
parameter estimator. Its general expression is [Ref.15]: 

^MUSIC^ p 
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where M is the number of complex sinusoids in the signal and 
v' are the eigenvectors that span the noise subspace 
(associated with the smaller eigenvalues). 

The discrete signal components will be orthogonal to 
any linear combination of the eigenvectors that form the noise 
base, and so we could expect singularities at the signal 
frequencies. In reality, we will not have pure singularities 
but the spectral estimates will peak. For our purpose of 
estimating a single frequency component, the procedure used 
was based on finding the roots of the averaged noise space 
eigenvectors (taken as polynomials of the complex 
exponential), whose phases correspond to the peak frequencies. 
The eigenvector associated with the root closest to the unit 
circle is the single frequency estimator. 

It can be proved that when we use an order 2 matrix in 
a single complex frequency component signal, the frequency 
estimator is equivalent to the Pulse Pair algorithm, and so 
the frequency estimate can be obtained directly from the 
autocorrelation function. 

The autocorrelation function for a single complex 
component signal superimposed in white noise will have the 
form: 

r^ik) =P. exp (j2nf-^ ) +p„6 (k) 
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and the corresponding 2"'^ order autocorrelation matrix will be: 




The eigenvector associated with the smaller eigenvalue (-p„) 
will be: 


i?(0) /2l^ 


\RiO) 
( 1 ) 


^2 


exp ( -j2E-^) 

^ S 


Therefore, the determination of the frequency estimate in this 
case will correspond to the determination of the phase 
associated with the first lag autocorrelation function 
estimate, as in the Pulse Pair estimator. 

This analytical solution of the 2"'^ order eigenvalue 
problem was implemented for performance comparisons, and a 
singular value decomposition was used for higher orders. The 
Modified Covariance Method was used to estimate the 
autocorrelation matrix at any order. 

B. MONTE-CARLO SIMULATION 

To evaluate the performance of these spectral estimators 
for the estimation of the frequency of a complex sinusoid 
superimposed in complex independent gaussian white noise, a 
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MONTE CARLO simulation was built within MATLAB environment. 
This simulation is based on 50 independent trials where the 
BIAS, VARIANCE and a TIME INDEX are evaluated at each step. 
The model has the input variables: 

1. Signal to Noise Ratio (scalar or vector) 

2. Number of Points (scalar or vector) 

3. Signal Frequency (scalar or vector) 

4. Type of Spectral Estimator 

The model output variables are: 

1. Bias 

2. Variance 

3. Time Index - proportional to the computation time 
and are stored at each step. 

The algorithm of the simulation program is as follows: 

1.SELECT TECHNIQUE 

1. Peak FFT 

2. Pulse Pair 

3. Zero Cross 

4. Modified Prony's Technique 

5. AR Model (burg algorithm) 

6. MV (capon's method) 

7. Music (eigenvalue decomposition'! 

2.SELECT PARAMETERS fincrements.initial and final values) 

1. Number points <-> SNR 

2. Frequency <-> SNR 
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3.Frequency <-> Number Points 


3.RUN MODEL 

1. Initialize variables 

trials=50 

test frequency=20 Hz 
sampling frequency=100 Hz 

2. Main loop 

generate noise 

generate new frequency 

generate signal 

estimate frequency 

compute estimate bias and variance 

compute processing time 

The complex independent gaussian white noise is generated 
by a function that uses different seeds for the real and 
imaginary parts at each step, in order to generate independent 
complex gaussian noise. Each noise component is then averaged 
from 5 different trials to obtain a pseudo gaussian noise. 
This number of 5 trials represents a compromise in the 
computation speed. The generated noise is normalized to 
guarantee a variance of 1. 

Model runs using a 20 Hz signal, a sampling frequency of 
100 Hz, over a -20 to 10 SNR and 16 to 128 points time series, 
were made for the following spectral estimators: 

1. Peak FFT 
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2. Pulse Pair 


3. Zero Cross 

4. AR Model order 4 

5. MUSIC order 4 

6. MUSIC order 3 

7. MUSIC order 2 

8. MUSIC order 2/Pulse Pair (fast algorithm) 

The selected plots of the bias and the logarithms of standard 
deviation are shown on the Figures 3.B.1-6. 

1. Discussion of Some Results 

The performance of an estimate can be assumed good 
when its bias is less then 1 Hz and the logic of the standard 
deviation is less then 0. Using this criteria the following 
tables show the SNR at which the bias of the several 
estimators is equal to 1 Hz, and the log standard deviation 
equal to 0: 

TABLE 1 - SNR at BIAS = 1 Hz 


Estimator 

20 Points 

40 Points 

60 Points 

80 Points 

Peak 

-4 

-5 

-7 

-8 

FFT 





Pulse 

-3 

-8 

-8 

-8 

Pair 
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(order 4) 
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(order 4) 
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Peak 

FFT 


Pulse 

Pair 
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Estimator 

20 Points 

40 Points 

60 Points 

80 Points 

MUSIC 

(order 4) 

10 

5 

2 

2 

MUSIC 

(order 2) 

5 

3 

2 

0 


2. Conclusions 

For all the estimators, after a sequence length (NP) 
between 20 and 64, the estimator variance becomes dependent 
only on the SNR. As a general rule except, for AR modeling, 
the variance is more restrictive then the bias. This means 
that we can get some improvement in the performance of this 
estimators if we take averaged measurements of the estimated 
frequency. The estimator that seems to work better for low SNR 
is the PEAK FFT followed by AR modeling (order 4), but AR 
modeling has poor processing time indexes (too slow). The 
MUSIC estimator is very sensitive to the noise contamination 
(SNR). The equivalent PULSE PAIR method and the MUSIC (order 
2 fast algorithm) do not show any significant improvement by 
changing the AUTOCORRELATION estimators, at least between the 
two that were used in this program (biased and the modified 
covariance estimators [Ref.l5]). 

The fastest algorithm is the MUSIC (order 2 fast 
algorithm) at all sequences lengths followed by the PULSE PAIR 
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riqnre 3.B.2 - Monte-Carlo Simulation: Pulse Pair algorithm, 
Bias and Variance plots. 
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Fiqiirc 3.B.5 - Monte-Carlo Simulation: MUSIC order 4 

algorithm, Bias and Variance plots 
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and ZERO CROSS for small sequences and the PEAK EFT for large 
sequences. 

The lowest SNR for which we can get unbiased estimates 
with a good variance is -10 when using a PEAK FFT with 128 
point sequences. Unfortunately this performance indicator does 
not show the discretezation error, which becomes 
overwhelmingly significant for small FFT lengths. We can try 
to reduce this error by using interpolation schemes, however 
this procedure will make this algorithm slower for a small 
number of points and is very susceptible to low SNR. For SNR 
between -5 and 0 we should use sequences of at least 32 to 64 
data points (PEAK FFT also), but these estimates will have 
severe discretization errors. For SNR between 0 and 5 we can 
use either the PEAK FFT or the PULSE PAIR/MUSIC 2 with 128 to 
32 sequence length. For SNR greater then 5 the best 
performance is given by the MUSIC (order 2 fast algorithm) or 
PULSE PAIR. If we want to emphasize the speed of the 
estimators, the best results are the ones given by the PULSE 
PAIR/MUSIC 2, requiring high SNR. 

Because the CDV package bistatic doppler frequency 
estimation is based on 32 points records, and the algorithms 
run in real time using C compiled code in a DSP, the most 
suitable algorithms selected from these simulations are the 
PEAK FFT for poor SNR, and the PULSE PAIR/MUSIC 2 algorithms 
for higher SNR. When considering discretization noise, the 
first choice for a spectral estimator to be used in the CDV 
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real time processing will be the PULSE PAIR estimator, because 
of its performance, simplicity and velocity. 
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IV. THE CDV MICROSTRUCTURE PACKAGE 


A. RECENT METHODS FOR TURBULENCE OBSERVATIONS 
1. Rapid Sampling Vertical Profiler (RSVP) 

The RSVP was developed by Caldwell, Dillon and Mourn in 
1987 [Ref.19] for sampling near microstructure parameters from 
a moving vessel permitting the study of the horizontal 
variability of upper layer profiles. The instrument makes 
near-microscale measurements (3 cm resolution) of temperature, 
salinity and turbulence dissipation down to 00 m when 
deployed from a ship moving up to 6 Knots. 

The profiler consists of a small retrievable probe at 
the end of a cable that is also a data link. Damping flaps are 
used to ensure uniform vertical velocity and orientation. The 
airfoil shear probe, used to estimate the TKE dissipation 
rate, is very sensitive to vibrational noise of the profiler, 
which may dominate the signal in many situations. The process 
by which it senses shear is based on the deflection of a piezo 
ceramic beam in a rubber jacket by cross-stream velocity. By 
careful instrument design, virtually all current measurements 
of ocean turbulent velocity have been made with this type of 
sensor. 

Mourn et.al. [Ref.10] performed upper ocean 
observations with this instrument where they were able to 
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characterize the fine structure of transition zone coastal 
eddys, and Dewey et.al. 1990 [Ref. 7] used this type of 
measurements to study the enhancement of ocean fronts by 
turbulent mixing. The interesting results obtained in both 
experiments show the importance and applicability of the RSVP 
in upper ocean observations- 

2. Pulse to Pulse Coherent Sonar 

A pulse to pulse coherent sonar was developed by 
Lohrmann et.al. 1990 [Ref.20], called the High Resolution 
Current Profiler. This instrument uses four 307 KHz 
transceivers, mounted 30“ off the vertical in 90° azimuthal 
intervals. Fifty range cells, covering a profiling distance of 
10 meters, measure radial current components along the 
transceiver beams. These measurements do not give the true 
velocity vector, but assuming a stationary and homogeneous 
turbulent field, allows an estimation of the Reynolds stresses 
and turbulent kinetic energy by separating the mean and 
fluctuating parts of the radial speeds. In this way, this 
instrument remotely senses the turbulent stresses without 
inducing any flow distortions and provides moderate vertical 
resolution. 

3. Bistatic Acoustic Doppler Profiler 

This instrument, developed by T.P. Stanton (1990), 
measures the explicit velocity vectors at a very high temporal 
and spatial sampling rate along one direction. The system 
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consists of four, 300 KHz fan-beam receiver transducers, 
radially separated from a central narrow beam-width 
transmitter transducer, which rapidly emit short duration 
acoustic pulses. Range gated doppler shifts are evaluated at 
each receiver using dual interpulse period, coherent doppler 
estimation methods. 

Simultaneous receptions at the four transducers 
provide an overdetermined estimate of the three velocity 
components at range bins along the ensonified volume, 
typically every 10 cm for 10 ro ranges. This system allows 
concurrent measurement of reynolds stress components and fine- 
scale shear out to ranges of 30 m. 

Velocity ambiguities from the coherent velocity 
estimation limit the unambiguous velocity range and profiling 
range. Unfortunately, the large frame for mounting the 
acoustic components restricts applications of the system. 

B. SPECIFICATIONS OF THE CDV SYSTEM 

The CDV Microstructure Package is a system developed by 
T.P. Stanton (1991) to measure ocean turbulence variables at 
a sub-centimeter resolution. It has 3 convergent acoustic 
transceiver beams operating at 5.296 MHz in a backscatter 
coherent mode, using both bistatic (0.6 cm focal cube) and 
monostatic operating modes. In the bistatic mode, one 
transducer is transmitting, while the other two are receiving 
the scattered energy from a focal volume. In the monostatic 
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mode the same transducer that transmits the signal receives 
the backscattered energy at range bins along the acoustic 
beam. 

The acoustic measurements provide the velocity vector, 
shear and kinetic energy estimation at a 44 Hz rate, while 
simultaneous microstructure measurements of temperature, 
conductivity, pressure, shear and the 3 components of the 
instrument acceleration are made. 

This package operates in a profiling mode at approximately 
25 cm/s descending or ascending rate, or fixed to a tripod 
facing down for bottom layer observations, or facing up for 
shallow water, near surface, observations. The package has 
also been attached to an ROV for semi-lagrangian upper ocean 
turbulence measurements. 

The digital data is transmitted to the surface by a 1.4 
Mbaud fiber optic data stream that is stored on a PCM/VCR 
system. The data stream is also passed to an AT&T digital 
signal processor (DSP) in a 386 based HP VECTRA for real time 
processing and graphical monitoring of selected channels. This 
pre-processed data is then stored on optical disks for further 
processing. A summary of the profiler package technical 
specifications follows: 
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LENGTH 

DIAMETER 

WEIGHT 

1.2 m 

12 cm 

7 Kg (in air) 


TABLE 4 - ACOUSTIC CHANNELS 


FREQUENCY 

POWER 

SAMPLING 

OUTPUT 



FREQUENCY 

FREQUENCY 

5.296 MHz 

5 Watt peak, in 

11.1 KHz (after 

44 Hz 


a 2 ° beamwidth 

demodulation) 



TABLE 5 - ANALOG CHANNELS 


SAMPLING 

TEMPERATURE 

CONDUCTIVITY 

OUTPUT 

FREQUENCY 

PRE- EMPHASIS 

PRE-EMPHASIS 

FREQUENCY 

1.4 KHZ 

0.339 Hz 

0.100 Hz 

44 Hz 

5.6 KHz (shear) 





C. SIGNAL PROCESSING OVERVIEW 

The acoustic inodes are controlled by a programmable 
sequencer that rotates the transmitting and receiving modes of 
the three transceivers, such that 3 monostatic and 3 bistatic 
modes have been selected every 20 ms. This acoustic data is 
received, complex demodulated to baseband and the complex (in 


64 



































phase and quadrature) channels sampled at 11.1 KHz, for the 
bistatic data and 88.2 KHz for each monostatic channel. 

The temperature and conductivity variables are pre¬ 
emphasized, in order to obtain an improved SNR at small scale 
and directly estimate gradient quantities, and sampled at 1.4 
KHz. The 3 acceleration channels and pressure are sampled at 
the same rate, while the shear pre-emphasized data is sampled 
at 5.6 KHz. 

The data stream comprises a set of 8 word records, such 
that the first two are used by the non-acoustic channels and 
the other 3 pairs correspond to the complex demodulated 
acoustic data. These 8 word records are contained in a 
sequence frame of 2051 words (corresponding to 22.7 ms of 
data) which are passed by DMA to the DSP for processing and to 
a PCM/VCR system for stor<ige. 

Once the DSP acquires the data, the data stream is 
interpreted and the real time processing is performed. 
Further processing, such as: 

1 - Water salinity computation 

2 - Water density computation 

3 - Velocity vector (u, v, w components) estimation 

4 - Spectral analysis of several variables 

5 - Reynolds stresses estimation 

6 - Scalar fluxes estimation 

7 - Dissipation estimation 

8 - T.K.E. estimation 
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is performed by external programs accessing this preprocessed 
data on the optical disks. The overall flow graph is shown in 
the figure 4.C.I. 

D. ACOUSTIC CHANNELS - VELOCITY ESTIMATIONS 

The velocity component estimation in the acoustic bistatic 
mode is based on the mean doppler shift of the backscattered 
acoustic energy resulting from the relative velocity of 
scatterers present in the focal volume. The real-time pre- 
processed acoustic data stored on the optical drive comprise 
the mean doppler shifts for the 6 different acoustic paths 
between pairs of transmitters and receivers. In the bistatic 
mode, the single transmitter/two receiver configuration is 
sequenced around the three transducers. In order to obtain the 
velocity components from this spooled data, the following 
steps are performed: 

1 - Clear data of large spikes using a reciprocal 
channel comparison and a 4*'’ order linear predictor, based on 
the previous 20 point autocorrelation function estimate. 

2 - Compute velocities from the doppler frequencies. 

3 - Transform estimates to an orthogonal frame of reference 
fixed to the instrument. 

4 - Rotate estimate to an absolute vertical reference frame. 

1. Doppler Velocities Estimation 

For two generic transducers, where transducer 1 is 
transmitting and transducer 2 is receiving, see figure 
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CDV Signal Processing overview. 

























4.D.1.1/ and assuming that the equipment is fixed to the 
bottom, we can derive the following doppler relationships 
[Ref.21]. If a pulse is transmitted at time t=0, it will 
arrive at BS at time t=ti, that for a sound velocity "c" will 
be given by: 

C-Vbs 


Now, if a second pulse is transmitted at time tj, it will 
arrive at BS at time t, and 



Figure 4.D.1.1 - Scheme of a Bistatic 2 transducer system. 
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Assuming that during the interval (Ojtj) we are transmitting 
at a frequency there will be wave fronts arriving at 
BS, during the interval t^-t^. Therefore, the frequency 
detected at BS will be: 


^BS ' ^ ^ 


= f ^ =»f = 

tx‘ -2 ^BS 


c-v, 


T1 

BS 


. f 


Cx 


This frequency corresponds to the frequency of the scattered 
signal at BS that will be send towards T2. Using the same 
procedure we arrive at the following expression for the 
frequency at T2: 


^ . f . f 

C+VfV C-*-Vps 


If we suppose that the velocity at BS has one component V^, 
see figure 4.D.1.2, then and f-n = ft^* Therefore 
the only velocity component producing a doppler shift will be 
along the direction defined by the axis bisection of the two 
beams i.e. V^. This component will have equal projections 
along the two beams, so that: 

^Bs = '^Bs 


Therefore 
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Figure 4.D.1.2 - Velocity projections scheme for a 2 

transducers bistatic system. 


I’'' =r. 
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- Af’■ 
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2 - • 


and 


V 



r . 


AX'- 


Because of the symmetric geometry of this system, this 
expression will also be true for the case when T2 is 
transmitting and T1 receiving. This fact creates redundancy in 
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Figure 4.D.2.1 - CDV Package, Orthogonal Reference Frame 

fixed to the instrument. 


the observed doppleL frequency, which provides increased 
confidence of the estimates by averaging, and allows the 
detection of single spurious wild points. Therefore the six 
acoustic transceiver modes, will measure three slanted 
velocity projections along the axes bisecting the acoustic 
beams. 

If the instrument was in motion we would have to add 
the instrument velocity components along these axis, in order 
to get absolute values. These operations are implemented by 
programs external to the real time processing. 






2. Velocities on an Orthogonal Reference Frame 

The Doppler velocity components must be transformed 
into an orthogonal frame before the Reynolds stresses or the 
scalar fluxes can be estimated. We can define the instrumental 
orthogonal reference frame as in figure 4.D.2.I. 

The scatterers velocity will be defined as: 

y .=Uji+Vjj*Wjk 


and the observed velocities will be: 




cosy 

0 

siny 




cosycosa 

cosycos (a--^) 

siny 

X 

' ,1- 


cosycosa 

-cosycos ) 

siny 



i . e. 


Therefore, this system may be solved in terms of the unknowns 
u, v and w by inverting the non-singular matrix "T". 

The next step is to rotate this velocity in the x-axis 
and y-axis in order to obtain a vertical z-axis. This can be 
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accomplished in several ways depending on each particular 
deployment. For the case when we have the system fixed to a 
tripod, we can look for the rotations that minimize the mean 
vertical velocity, or just measure the platform orientation 
using the 3 axis accelerometer as tiltmeters. 

E. NON-ACOUSTIC DATA 

The non-acoustic channels on the CDV package are: 3 
component acceleration, pressure, temperature,conductivity and 
shear. The real time processing for these channels was coded 
in C language and compiled on the DSP together with the 
Assembler routines that access the CDV data stream in memory. 
This DSP program is controlled by another C routine running on 
the Host computer, that interacts with the user, reads the 
proper setup files, receives the data from the DSP, plots 
selected variables of the pre-processed data and sends the 
output to an optical mass storage unit. The details on the 
processing of each particular channel are described next. 

TEMPERATURE AND CONDUCTIVITY 

The temperature variable (T) is measured using a high 
frequency response FP 07 thermistor, and conductivity (C) 
using a 2 electrode AC conductivity probe. The analog 
temperature and conductivity are pre-emphasized at 20 dB per 
decade for frequencies greater than 0.339 Hz and 0.100 Hz 
respectively, for SNR optimization. They are sampled at 1.4 
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KHz and the multiplexed data stream is sent to the DSP by DMA 
assembly input routine. 

A DSP function demultiplexes variables from the data 
frame to separate selected variables from the stream. A first 
difference criteria with a threshold s=t by a setup file, is 
used to detect spikes on the data. Single bad points are 
substiif.uted by the previous valid point. Deterministic 
coherent noise over each frame period (32 words) is subtracted 
from each data cycle. This coherent noise is computed from the 
mean of approximately 5000 frames and is entered from an input 
file. The raw T and C gradient data are then squared and 
accumulated on 1 second rotary buffers for estimation of 
gradient variances. Both the T and C data are also de- 
emphasized by single pole digital Butterworth filters and low 
pass (LP) filtered by a second order Butterworth filter with 
a cutoff frequency at 20 Hz. These operations are done in 
sequence by different functions to avoid round off errors. 
Finally the de-emphasized data is decimated to 44 Hz and the 
variances are computed based on the mean of the 1 second 
rotary buffers. These 44 Hz data are sent every second to a 
HOST program running on the HP Vectra. 




Acceleration is observed using 3 orthogonal solid 
state DC response accelerometers in the head of the CDV. 
Pressure is sensed using a strain gauge pressure sensor. Both 
Acceleration and Pressure are processed in a similar manner. 









The DSP function that demultiplexes the input buffer, also 
separates these 4 channels from the data stream. The same 
deglitching procedure is used concurrently with a threshold of 
200 raw counts for the four channels. There is negligible 
observed coherent noise on these channels. The data is then LP 
filtered by a third order Butterworth with a cutoff frequency 
at 20 Hz, decimated to 44 Hz and sent to the HOST program 
every 1 second for monitoring and storage. 

SHEAR 

To measure the velocity shear, cross-velocity 
variations are measured with a high sensitivity airfoil probe. 
Because of its nature, this sensor is easily noise 
contaminated by variations in the instrument velocity or 
spurious vibrations. Also, it only works when the instrument 
is in motion and senses only velocity changes perpendicular to 
the mean sensor velocity. 

The shear probe is sampled at 5.6 KHz and the data 
send to DSP by DMA assembly input routine. The same DSP 
function separates the shear from the data stream and applies 
a 200 count, first difference deglitching criteria. Any 
deterministic coherent noise is subtracted from each data 
cycle. Next, the pre-emphasized (gradient) data is squared and 
accumulated in 1 second rotary buffers, and the variance is 
computed at a 44 Hz rate and sent to the HOST program. 
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F. CALIBRATION AND ERROR ANALYSIS 


The temperature data was previously calibrated using a 
third order polynomial fit of the form: 


TEMP=A. COUNTS^ *B. COUNTS'^ *C. COUNTS 
where 
/9=10.662337 
fi=-3.957935,10'^ 

C=l.216833.10-^ 
£>=-1.345511. 


where TEMP represents temperature values in Degrees Celsius 
and COUNTS the CDV outputed values. 

The pressure sensor was calibrated by recording the 
pressure values obtained from water columns from 0.10 meter to 
2.00 meter, 0.10 meter apart. The calibration transformation 
was then computed using a least squares fit to a straight line 
by solving a system of 20 equations on 2 unknowns. The 
resulting calibration is: 

DEPTH=0 .000191964. COUNTS-0 .0218414 

where DEPTH is the water column height in meters and COUNTS 
the pressure channel output values. The mean square error is 
0.0012 m. The calibration was performed with pure water, so 
that in order to compute the depth in salty water, this values 
must be multiplied by a factor related to the density anomaly. 
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that for practical use can be 0.974 (assuming a density of 
1026 Kg/m") . 

The acceleration channels were calibrated using the fact 
that when the equipment is fixed the mean observed value 
should be the gravitational acceleration. The package was 
fixed and rotated from -20° to +20°, in 5° intervals on two 
quasi-orthogonal planes, with one of the planes containing one 
of the horizontal accelerometers. Least square fits to 
straight lines were performed, resulting in the following 
calibration: 

ACC^= 0.00108549 . . 0324 

ACC^= 0.00119582 . CO£WTSy-19.967 3 
ACC-=-0 .00110622 . COa2VT5-+16.9346 


where ACCi is the acceleration in m/s^ and COUNTSj the outputed 
quantities. The observed mean square errors are 0.0532 m/s^, 
0.0210 m/s^ and 0.0634 m/s^ in the X, Y and Z channels 
respectively. 

To evaluate the performance of the equipment, two 5 
minutes data sets were selected from the WHARF #2 Experiment 
at the August 2 1991. This experiment will be described in 
more detail in the next chapter. The calibrated temperature 
values are very closed to those observed using a mercury 
thermometer bucket temperature. The accelerometer mean values 
are reasonably close to the expected. Because there was no in 


77 







situ calibration, only mean quantities can be analyzed. A more 
detailed description of these data will be presented in the 
next chapter. This experiment was performed with the 
instrument fixed to the bottom, so no shear probe data was 
analyzed. No pressure data was available from this deployment. 

The first step in the analysis is to look at the pre- 
processed doppler data. A sample of two bistatic reciprocal 
channels (frequency, power and bandwidth), and calibrated 
temperature data is shown in the figures 4.F.1-2. 

A careful analysis was required for the velocity 
estimation, because there was no previous experience with this 
method. It was necessary to validate all the processing steps 
(both electronic and software). 

The reciprocal channel doppler frequency estimates are 
generally similar, see figure 4.F.1, except for individual 
spikes, or in regions of very low signal to noise ratio (SNR) 
which are also seen by a large bandwidth. Some of these 
discrepancies will be attenuated by using a variable pulse 
repetition interval in future deployments, eliminating 
coherent additive interference. 

To increase the confidence and deglitch the acoustic data, 
an adaptative procedure has been developed. The probability 
distribution function (PDF) of the differences histogram (see 
in the Figure 4.F.3), shows that the differences between two 
consecutive values are smaller than 70 Hz (0.01 m/s) 75% of 
the time. Also, the scatter plots between the first 
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Figure 4.F.1 - Reciprocal Bistatic channels, dopipler 

frequencies, received powers and bandwidth. 
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Figure 4.F.2 - Calibrated Temperature data. Simultaneous to 
the data shown in figure 4.F.1 . 

differences and the received powers (Figure 4.F.4) shows that 
the bigger differences are mainly associated with low received 
powers, corresponding to the cases when noise from additive 
interference and system noise may be significant. We have two 
independent realizations of the same random variable 
(reciprocal channels) that can be used either to increase the 
confidence of the estimates or to validate them. 

To deglitch the data and obtain a clean averaged time 
scries of the doppler frequencies, an intelligent adaptative 
algorithm must be used based on the surrounding points. A 
procedure was implemented to compare the signal at each time 
with a predicted value obtained from a 4’^ order linear 
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Figure 4.F.3 - Reciprocal Bistatic channels, first 

differences histogram. 


predictor based on the previous 20 point autocorrelation 
fvjnction. Whenever the difference was greater than four times 
the standard deviation (also estimated from the previous 20 
points), that point was assumed non-valid. When both channels 
were assumed bad, or if a first difference fixed criteria was 
violated, the predicted value was used. We can expect an 
improvement in the SNR of the estimates by using this method. 
.Some examples of how this procedure works, are shown in the 
figure 4.F.5. 

To transform these deglitched doppler shifts to veloci ty 
components in an orthogonal frame fixed to the instrument, a 
linear transformation must be performed. This transformation 


«1 










between first differences and received power. 

will not distribute the velocity variance evenly through the 
3 orthogonal components of velocity, giving more weight to the 
component along the instrument vertical; this componeiit will 
also have a greater statistical confidence. 
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DEGLITCHING FUNCTION - thresholds and errors 
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Figure 4.F.5 - Bistatic channels deglitching function. 
Whenever the difference curves go higher than the thre.shold 
curves, in the lower plot, that point was assumed as bad data. 
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V. •nJRBULENCE OBSERVATIONS USING THE CDV PACKAGE 

A. WHARF 12 EXPERIMENT 

The purpose of this experiment was to asses the 
performance of the CDV package in oceanic backscatter 
conditions, test its ability to measure ocean fluxes and 
surface wave fields, and to study the forcing mechanisms in 
ocean surface layers. 

The CDV package was deployed off the Monterey Wharf #2 on 
2 and 3 August 1991. The tripod-mounted CDV was placed on the 
sand-bed with the top of the instrument at depths ranging from 
0 to 0.75 meter. The depth changed with the tide, and the 
total water depths ranged between 1.7 to 2.5 meter. 

The package was controlled and monitored from a van parked 
on the wharf, from which observations of air temperature 
(moist and dry), bucket sea surface temperature, short and 
long wave radiation and wind direction and speed. The 
equipment deployment is sketched on figure 5.A.I. 

The air temperature was measured using a mercury 
psychrometer and the sea surface temperature using a protected 
mercury thermometer. The wind speed and direction was recorded 
using a MRI anemometer and wind vane, while the short and long 
wave radiation was obtained from an Epply radiometer. 
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Figure 5.A.1 - Wharf #2 Experiment set-up 








The environmental conditions were characterized by an 
overcast sky, moderate to weak winds from WNW. Small amplitude 
swell from the North was being reflected from the solid walls 
below the wharf. Due to the barrier created by the wharf walls 
on the wind field, the measured winds did not correspond well 
to those over the CDV package. These light wind conditions 
were not ideal to evaluate the overall performance of the 
instrument, because we would not expect to find high 
turbulence levels in the high frequency range. Nevertheless, 
because we understand the physics of the linear wave fields 
and we do not have lateral boundary effects (as we would 
encounter in a tank experiment), we were able to measure the 
performance of the package at the low and medium frequency 
range using the wave field signal, and to estimate the noise 
floor at higher frequencies. 

The CDV data were recorded on VCR tapes (labeled WCDVl to 
WCDV7), from 1100 until 2300 on 2 Aug 1991, covering the 
afternoon high tide at 1545 pm, and from 0200 to 0800 at 3 
Aug, covering the morning high tide at 0515. The short and 
long wave radiation and the wind data were recorded 
continuously from 1100 am of 2 Aug until 0800 am of 3 Aug, 
sampled at 8 Hz and spooled by an HP 3 32 computer to an 
optical drive for further processing. The air and bucket 
temperature were recorded every hour from 1600 2 Aug to 0800 
3 Aug. 
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After a careful analysis, only the afternoon CDV data from 
1300 to 1800 show high data quality, while the morning data 
show large data gaps due to the fact that the equipment was 
only partially submerged. 

Following DSP preprocessing, all the data were spooled to 
an HP VECTRA 386 and the CDV bistatic data to a SUN SPARC 
STATION for further processing. 

After calibration, the 8 Hz sampled wind and radiation 
data were 10 min averaged and decimated at a 1/min sampling 
rate. The CDV sea temperature was calibrated and used for 
cross-spectral analysis and computation of fluxes over 20 
minutes averages. The bistatic doppler frequency estimates 
were deglitched using the adaptative procedure described in 
the previous chapter, and used to estimate the velocity vector 
in an instrumental reference frame. The instrumental 
velocities were then properly rotated to a vertical reference 
frame. All the processed data (intermediate and final), 
spectral and flux estimates and output plots were stored on an 
optical disk. 

B. METEOROLOGICAL DATA 

This section describes and analyzes the tide, wind, 
radiation, air and sea surface temperature data. It will focus 
only on the period from 1100 to 2000 of 2 Aug 1991, 
simultaneous with the useable CDV data. 
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As is shown in the Figure 5.B.1, within this period there 
was a high tide at 1545 pm. The tide height determines the 
depth of the CDV measured quantities; so in Figure 5.B.1, the 
expected mean depth of the top of the CDV is shown, assuming 
that at 1100 it was at a 0 m depth. Also shown as vertical 
straight lines are the times of the evaluated ocean fluxes. 

Due to problems in the radiation sensors, long wave 
radiation data was not available. This fact restricts our 
ability to compute the surface buoyancy fluxes; however, we 
can use the short wave radiation and wind stress as indicators 



CDV package at 2 August 1991. 


88 





of the forcing mechanisms in the ocean upper layer. The wind 
data does not accurately reflect wind stress at the CDV 
location due to the barrier created by the wharf walls, but it 
will show the maximum expected mechanical forcing of the ocean 
surface. The time evolution of these variables is shown in 
Figure 5.B.2, for the period from 1200 to 2000 pm at 2 Aug 
1991, based on measurements made from the top of the wharf. 

These plots show that the wind stress was decreasing 
slightly after 1500. Due to the overcast conditions during 
almost all the period of the experiment, we should not expect 
a large variability on the long wave radiation, so the changes 
in the surface buoyancy flux should be determined primarily by 
the short wave radiation and by the surface layer stability. 
The short wave radiation followed a cloud covered diurnal 
cycle increasing until 1500, and decreased thereafter. The 
peak shown at approximately 1420 corresponds to a short period 
of time of partially clear sky. From these conditions, we may 
expect a decrease in the weak mechanical steering due to 
decreasing winds and a reduction in the stability of the upper 
layer and due to a reduction of solar heating. Therefore, we 
should not expect to find high levels of turbulent energy in 
the form of scalar fluxes and Reynolds stresses, i.e. we 
should not expect to observe an active mixing layer, so that 
the ocean conditions should be dominated by surface gravity 
wave fields. 
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Figure 5.B.2 - Wharf #2 Experiment: Wind stress and short 
wave radiation, from 1200 to 2000 at 2 Aug 1991. 


C. CDV DATA 

The data between 1500 and 1800 of 2 Aug 1991 were analyzed 
in more detail and the ocean fluxes computed using an eddy 
flux correlation method [Ref.2]. The first step was to select 
20 minutes ensembles every 1 hour. These data sets were then 
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spooled to an optical drive, and the data analyzed in more 
detail using C programs and MATLAB software. 

After conversion to engineering units, the acceleration 
and temperature data were averaged over all the 20 minute 
ensembles; the tilt angles were computed and the mean 
temperatures compared with the observed bucket temperature. 
TABLE 6 - CDV AVERAGED DATA 


Hour (2 Aug) 

15:30 

16:30 

17:30 

18:30 

Acceleration-X 

(m/s^) 

0.1626 

0.1588 

0.1630 

0.1553 

Acceleration-y 

(m/s^) 

-0.6645 

-0.6647 

-0.6718 

-0.6787 

Acceleration-z 

(m/s") 

-9.8278 

-9.8330 

-9.8277 

-9.8288 

Direction-Deg 

(ref.to x-axis) 

-76.25° 

-76.56° 

-76.36° 

-77.10° 

Tilt (Deg) 

(from vertical) 

3.98° 

3.976° 

4.02° 

4.05° 

Mean Temperature 

(bucket) 

17.1°C 

16.8°C 

16.3°C 

16.0°C 

Mean Temperature 

(CDV sensor) 

16.24°C 

16.23°C 

16.17°C 

15.82°C 
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Hour (2 Aug) 

15:30 

16:30 

17:30 

18:30 

CDV Depth 

0.74m 

0.71m 

0.64m 

0.36m 

(Tide Table) 






We can see that over the full deployment there was a 
constant tilt of 4 degrees from the vertical, and an offset of 
-76.6 degrees from the x-axis, defined as the bisection of 
transducers 1 and 2. These values were used to rotate the 
velocity data to obtain a vertical reference frame. Examples 
of the measured velocities are shown in Figure 5.C.I. The 
small magnitude and variability of these values is consistent 
with the observed conditions of low or nonexistent turbulence 
and small amplitude wave field. 

The mean temperature compares well with the bucket 
temperature, but the 44 Hz measured temperature, shows some 
high gradients (as is shown in Figure 5.C.2) not correlated 
with the vertical velocities, suggesting high horizontal 
stratification. The strong and sharp variations of 
temperature, by itself, suggest horizontal advection of warm 
patches of water. This hypothesis is further discussed in the 
next section, when the temperature fluxes are discussed. 

D. OCEAN FLUXES AND STRATIFICATION 

Figures 5.D.1 to 4, show the spectral evolution of 
temperature and of the 3 velocity components. The frequency 
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density units of all spectral quantities are such that the 
integration over the full frequency range gives the total 
variance. These spectra are characterized by a broad peak, 
more evident in the velocity spectra, at frequencies between 
0.2 and 0.6 Hz. They represent short period, surface gravity 
waves, with wavelengths of the order of 6 m with e-folding 



at 44 Hz. 

depths of the order of 1 m (i.e. deep water waves). This fact 
is consistent with the observed similar energy content in the 
vertical and horizontal velocity components (10‘* {m/s)*/Hz). 
The broad band shape of the spectrum suggests that this energy 
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is due to wind waves starting to shoal, probably generated in 
the near field. 

There are also three peaks, possibly associated with low 
frequency swell or seiche, at the frequencies of: 2.5x10'^ HZ, 
4.5x10'' Hz and 7x10'^ Hz: these spectral peaks can be seen in 
the figure 5.D.5, which shows the high resolution spectra of 
the variables at 1530. These peaks are more noticeable in the 
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Figure 5.C.2 - Wharf ^2 Experiment: Temperature data. 


sampled at 44 Hz. 

spectra of the horizontal velocity, significant to a 80% 
confidence level, and are highly attenuated in the vertical 
velocity spectra as we can expect for a shallow water wave. 
These harmonically related frequencies suggest weak non-linear 
interactions (resonance), with energy being transferred, 
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Figure 5.D.3 - Wharf il2 Experiment - Spectral evolution of 
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Figure 5.D.5 - Wharf #2 Experiment - High resolution spectra 
at 15:30, 
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between the different waves components. If ressonance is 
occuring, these components will contribute to the presence of 
stresses in the upper layer. 

Besides telling us that the velocity field was dominated 
by the surface gravity waves, the low energy levels at high 
frequencies also shows the noise floor of the instrument with 
energy density on the order of 10“’ (m/s)^/Hz in the velocity 
spectra, and 10‘^° (Deg C)^/Hz in the temperature spectra. This 
value represents noise significant velocities of the order of 
0.005 m/s, that are smaller than typical velocities associated 
with active oceanic turbulence, typically of the order of 0.01 
m/s [Ref.22]. Therefore, for turbulent flux observations we 
should expect to obtain good SNR performances with this 
instrument. 

Integrating the velocity spectra, we obtain the kinetic 
energy for each ensemble. These values and their evolution 
over the 4 analyzed hours are important in the 
characterization of the observed fields. We can see that the 
total horizontal kinetic energy was significantly larger than 
the vertical, due to the nature of the wave fields. 
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HOURS 

15:30 

16:30 

17:30 

18:30 

u-component 

(mVs^) 

0.0044 

0.0049 

0.0062 

0.0119 

v-component 

(mVs") 

0.0107 

0.0118 

0.0149 

0.0238 


In order to compute and interpret the fluxes of 
temperature and the Reynolds stresses, the co-spectra, phase 
spectra and coherence were computed for the several 
combinations of variables. Using Parseval theorem [Ref.14], 
integrating the co-spectra over frequency provides estimates 
of the scalar fluxes or the Reynolds stresses. Therefore, the 
co-spectrum will tell us how these fluxes are distributed over 
frequency. In figures 5.D.6, 7.a, 7.b and 7.c, we can see the 
evolution of the co-spectrum between the perturbed quantities 
of the temperature (t') and vertical velocity component (w') 
hereafter designated by t'w', and between the three velocity 
components u'w', v'w' and u'v'. The low values in the co¬ 
spectrum and the negligible energy above 1 Hz are 
characteristic of a linear wave field with low turbulence 
levels. 

The negative peak found at the very low frequency range in 
the co-spectrum of the horizontal velocity components. Figure 
5.D.7.a, confirms the presence of non-linear interactions 
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Figure 5.D.6 - Wharf it? Experiment - Time evolution of the co¬ 
spectrum, between temperature and vertical velocity 
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Figure 5.D.7.C - Wharf #2 Experiment - Tire evolution of the 
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between harmonically related wave components [Ref.23], as 
suggested by the velocity spectra in Figure 5.D.5. 

The co-spectra of temperature and vertical component of 
the velocity, shows that the thermal energy flux is negative 
in the short period wave range (deep water waves), and 
positive over the long period waves range at the 1530 and 1630 
hours. Later energy starts cascading into the smaller scales, 
where previously the co-spectral density was almost zero 
(Figure 5.D.6). This is partially consistent with an increase 
in the horizontal momentum flux, though it remains restricted 
in the low frequency range. No significant changes are seen in 
the vertical momentum flux. This picture suggests that after 
the 1730, due to the shallow depth of the instrument and to 
the reduction of insolation generated stratification, 
turbulent levels tend to increase. 

After integrating equivalent estimates for each data set, 
the time evolution of the temperature fluxes and stresses were 
calculated: 


TABLE 8 - TEMPERATURE AND MOMENTUM FLUXES 


HOURS 

15:30 

16:30 

17:30 

18:30 

<u'v'> 

(roVs^) 

-0.0046 

-0.0051 

-0.0063 

-0.0108 

<u'w'> 

(mVs") 

-0.0001 

-0.0000 

0.0002 

0.0015 
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HOURS 

15:30 

16:30 

17:30 

18:30 

<v'w'> 

(mVs") 

0.0011 

0.0014 

0.0020 

0.0012 

<t'w^> 

-2.2277 

-1.0943 

-2.5964 

-10.5082 

(m.°C/s) 

XlO'^ 

XlO'^ 

XlO'^ 

Xl0‘* 


This table shows that the Reynolds stresses and temperature 
fluxes slightly increased with time, as stratification was 
decreasing. 

The phase spectra and coherence are tools that help the 
interpretation of the co-spectrum. The phase spectra quantify 
the phase lags between variables, showing departures from the 
in-phase components shown in the co-spectra. These spectra are 
useful in the identification of quadrature conditions, 
characterizing linear wave components, and in the analysis of 
non-linear characteristics. The coherence spectrum shows 
whether cross-spectral components have constant phase 
relationships in time. For linear wave fields, coherence 
should be high, near 1 (no fluxes or stresses), while for 
turbulent or non-linear fields the coherence should have low 
values. Examples of these estimates are shown in the figures 
5.D.8 and 9 for the temperature vertical flux and Reynolds 
stresses at 1530. The features shown in these figures are 
similar to the characteristics of the other data sets, with 
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velocity components and temperature at 15:30. 
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Figure 5.1).9 - Wharf il7. Experiment - Coherence of velocity 
components and temperature at 15:30. 


109 







quadrature phase and high coherence through the spectral peaks 
range, as expected for a linear wave field. 

The low correlation between temperature and vertical 
velocity hint that we cannot explain the observed temperature 
variability only from the vertical stratification, in which 
water is pumped up and down past the sensor following the wave 
orbital motion. On other hand, if we assume that horizontal 
stratification exists, then the horizontal velocities may 
explain this variability, with patches of warmer or colder 
water being advected by the mean horizontal flow and wave 
orbits. To test this hypothesis, the horizontal temperature 
fluxes were computed as is shown in the figure 5.D.10 for the 
1530 data set. In fact, they do not show a stronger 
correlation; however, the probable high variability on the 
horizontal density gradient will introduce an additional 
effect on these estimates. These estimates may also be 
misleading if not computed from sufficiently long averaging 
intervals as the sensor was stationary in the field. 

The coherence spectra( Figure 5.D.9) again reinforce the 
fact that the wave regime at the middle frequencies has a 
linear character, showing high coherence between v and w 
within this range. These characteristics are maintained 
throughout the observation period. Similarly, the coherence 
between the horizontal components of velocity is more 
significant in the low frequency wave range, but with slightly 
lower values. 
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Figure 5.D.10 - Wharf §2 Experiment - Horizontal fluxes of 
temperature at 15:30. 
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VI. CONCLUSIONS 


The need to explicitly observe turbulent processes arises 
from a wide class of problems in physical oceanography, as 
described in the review in chapter II. The CDV package was 
developed to simultaneously measure temperature, conductivity, 
pressure, shear and the 3 components of velocity at sub¬ 
centimeter scales. 

The unique measurement techniques used in this instrument 
require special data processing and analysis methods. The 
three component velocity estimations rely on a doppler 
frequency spectral estimator, working in real-time, that must 
use a small number of coherently sampled data points. In order 
to choose the meet suitable a_gorithm, and to assess its 
physical limitations, Monte-Carlo simulations were performed 
using several possible spectral estimators. From these 
simulations, the Pulse Pair algorithm, or equivalently the 
MUSIC order 2 algorithm, showed the best performance, followed 
by the Peak FFT algorithm. The FFT estimator, corresponds to 
a maximum likelihood estimator, but suffers from the fact that 
for a small number of points, the discretization errors are 
very dominant, and it's performance is seriously degraded. 
Using these results, the Pulse Pair estimator was implemented, 
together with the pre-processing of the other measured 
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quantities, on a Digital Signal Processor for real-time pre¬ 
processing and review. 

To test the CDV software and hardware, the equipment was 
deployed for the first time off Monterey Wharf #2, on the 2 
and 3 of August 1991, along with measured environmental 
conditions. This first analysis showed that the performance of 
the bistatic velocity estimation is below that expected, with 
some conspicuous spikes superimposed rn valid data. By 
introducing a change in the transmitting cycles of the 
transducers this problem should be overcome in future 
deployments. The calibration of the other scalar variables and 
acceleration gave good results. To deglitch the bistatic 
doppler frequencies and increase its confidence, we used the 
redundancy in the velocity components observations, with two 
independent estimates from each of the two reciprocal acoustic 
bistatic paths. An adaptative procedure was developed, relying 
on a one lag linear-predictor and predictively choosing the 
best estimate of the doppler frequency. This procedure 
significantly improved the SNR of the velocity estimates. 

The CDV Wharf data shows a surface gravity wave field with 
energy in two different ranges, superimposing near field 
generated high frequency wind waves with a low frequency 
swell. The wind waves sho«» the expected characteristics of 
deep water waves while the low frequency waves clearly have 
shallow water waves characteristics. The high resolution 
spectra and the co-spectra show that resonance interactions 


113 






may be occurring at the low frequency range. The high 
frequency wind waves show a stronger coherence and do not 
contribute significantly to the momentum flux, but 
nevertheless the major contribution to the temperature flux 
occur over these frequency ranges. 

The fact that the observed temperature variability did not 
correlate well with the measured vertical velocity, even at 
the deep water wave frequency range, suggest that an irregular 
horizontal stratification existed, with pools of warm or cold 
water being advected past the sensor, by the wave orbits and 
mean flow. 

From the vertical velocity and temperature co-spectra 
analysis, we also observed an increase in the temperature flux 
associated with a shallow depth and reduction in the vertical 
stratification, following the decrease in the solar heating 
after 1730. However, due to remaining significant 
stratification no significant turbulence develops, as we can 
observe in the co-spectra of the velocity components. 

At the high frequency ranges, we did not found significant 
energy due to the absence of strong turbulence, allowing the 
instrumentation noise floor to be seen. This noise floor 
represents velocities in the order of 0.005 m/s, values below 
the ones that we should expect to find during strong turbulent 
conditions. Because this data was the first that was analyzed 
with this equipment, we expect in future deployments to 
improve the SNR in the velocity measurements by an order of 
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magnitude by reducing the acoustic noise component 5 to 10 
times and altering acoustic sampling methods. However, the 
existing noise floor is already below the minimum requirements 
for studying energetic turbulence. 

The CDV data seems consistent with the observed 
environmental conditions, confirming the expected capabilities 
of the instrument to observe real data, including wind wave 
fields, and showing the noise limitations that we should 
expect when observing turbulence fields. However, in order to 
further study its precision, this data should also be compared 
with simultaneous observations of the ocean variables measured 
by independent instruments. The CDV package will be able to 
perform simultaneous measurements of oceanographic variables 
at the sub-centimeter scales, determining all terms in the 
Reynolds equations. This should improve our characterization 
of turbulent fields, with applications in a variety of 
oceanographic turbulence problems. These applications will be 
explored using the instrument in a loose-tethered profiling 
mode. 
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